########################################################################
################# Describes DIPCON network 1970-2010 ###################
########################################################################

########################################################################
# Install packages
detach(package:statnet, unload=T)
install.packages("igraph")
install.packages("plyr")
install.packages("xtable")
install.packages("ggplot2")
########################################################################

# set working directory
setwd("~/Dropbox/ISQ_ms/5 Final/Replication")

# set savedir
# SAVEDIR <- "~/Dropbox/ISQ_ms/5 Final/Replication"

# load libraries
library(igraph)
library(plyr)
library(xtable)
library(ggplot2)

# read in data
# read in adjacency matrices
dipcon1970 <- read.csv("DIPCON1970_3.0.csv",row.names=1) # 134 states
dipcon1975 <- read.csv("DIPCON1975_3.0.csv",row.names=1) # 148 states
dipcon1980 <- read.csv("DIPCON1980_3.0.csv",row.names=1) # 157 states
dipcon1985 <- read.csv("DIPCON1985_3.0.csv",row.names=1) # 162 states
dipcon1990 <- read.csv("DIPCON1990_3.0.csv",row.names=1) # 164 states
dipcon1995 <- read.csv("DIPCON1995_3.0.csv",row.names=1) # 183 states
dipcon2000 <- read.csv("DIPCON2000_3.0.csv",row.names=1) # 187 states
dipcon2005 <- read.csv("DIPCON2005_3.0.csv",row.names=1) # 188 states
dipcon2010 <- read.csv("DIPCON2010_3.0.csv",row.names=1) # 190 states

# create graphs
g1970 <- graph.adjacency(as.matrix(dipcon1970), mode="directed",weighted=NULL)
g1975 <- graph.adjacency(as.matrix(dipcon1975), mode="directed",weighted=NULL)
g1980 <- graph.adjacency(as.matrix(dipcon1980), mode="directed",weighted=NULL)
g1985 <- graph.adjacency(as.matrix(dipcon1985), mode="directed",weighted=NULL)
g1990 <- graph.adjacency(as.matrix(dipcon1990), mode="directed",weighted=NULL)
g1995 <- graph.adjacency(as.matrix(dipcon1995), mode="directed",weighted=NULL)
g2000 <- graph.adjacency(as.matrix(dipcon2000), mode="directed",weighted=NULL)
g2005 <- graph.adjacency(as.matrix(dipcon2005), mode="directed",weighted=NULL)
g2010 <- graph.adjacency(as.matrix(dipcon2010), mode="directed",weighted=NULL)

# check graphs
g1970
g1975
g1980
g1985
g1990
g1995
g2000
g2005
g2010


### Online Appendix

### Overall Network Structure

### Table A1. Network Statistics

# number of nodes
number_states <- c(vcount(g1970), vcount(g1975), vcount(g1980), vcount(g1985), vcount(g1990), vcount(g1995), 
                   vcount(g2000), vcount(g2005), vcount(g2010))

# number of isolates
number_isolates <- c(sum(degree(g1970)==0), sum(degree(g1975)==0), sum(degree(g1980)==0), 
                     sum(degree(g1985)==0), sum(degree(g1990)==0), sum(degree(g1995)==0), 
                     sum(degree(g2000)==0), sum(degree(g2005)==0), sum(degree(g2010)==0))

# number of ties
number_ties <- c(ecount(g1970), ecount(g1975), ecount(g1980), ecount(g1985), ecount(g1990), ecount(g1995), 
                 ecount(g2000), ecount(g2005), ecount(g2010))
# density
density <- c(graph.density(g1970), graph.density(g1975), graph.density(g1980), graph.density(g1985), 
             graph.density(g1990), graph.density(g1995), graph.density(g2000), graph.density(g2005))

# proportion of symmetric dyads
sym_dyads <- c(1-(dyad.census(g1970)$asym/(vcount(g1970)*(vcount(g1970)-1)/2)),
               1-(dyad.census(g1975)$asym/(vcount(g1975)*(vcount(g1975)-1)/2)),
               1-(dyad.census(g1980)$asym/(vcount(g1980)*(vcount(g1980)-1)/2)),
               1-(dyad.census(g1985)$asym/(vcount(g1985)*(vcount(g1985)-1)/2)),
               1-(dyad.census(g1990)$asym/(vcount(g1990)*(vcount(g1990)-1)/2)),
               1-(dyad.census(g1995)$asym/(vcount(g1995)*(vcount(g1995)-1)/2)),
               1-(dyad.census(g2000)$asym/(vcount(g2000)*(vcount(g2000)-1)/2)),
               1-(dyad.census(g2005)$asym/(vcount(g2005)*(vcount(g2005)-1)/2)),
               1-(dyad.census(g2010)$asym/(vcount(g2010)*(vcount(g2010)/2))))

# global clustering coefficient
clustering <- c(transitivity(g1970), transitivity(g1975), transitivity(g1980), transitivity(g1985), 
                transitivity(g1990), transitivity(g1995), transitivity(g2000), transitivity(g2005), 
                transitivity(g2010)) 


### Table A2. Embassies Opened and Closed Relative to the Previous Period (1975-2010)

# check how many embassies were created between 1995 and 2005

new75 <- difference(g1975, g1970)
lost75 <- difference(g1970, g1975)

new80 <- difference(g1980, g1975)
lost80 <- difference(g1975, g1980)

new85 <- difference(g1985, g1980)
lost85 <- difference(g1980, g1985)

new90 <- difference(g1990, g1985)
lost90 <- difference(g1985, g1990)

new95 <- difference(g1995, g1990)
lost95 <- difference(g1990, g1995)

new00 <- difference(g2000, g1995)
lost00 <- difference(g1995, g2000)

new05 <- difference(g2005, g2000)
lost05 <- difference(g2000, g2005)

new10 <- difference(g2010, g2005)
lost10 <- difference(g2005, g2010)

newlosttab <- data.frame(matrix(c(ecount(new75), ecount(new75)/ecount(g1970)*100, 
                                  ecount(lost75), ecount(lost75)/ecount(g1970)*100,
                                  ecount(new80), ecount(new80)/ecount(g1975)*100,
                                  ecount(lost80), ecount(lost80)/ecount(g1975)*100,
                                  ecount(new85), ecount(new85)/ecount(g1980)*100,
                                  ecount(lost85), ecount(lost85)/ecount(g1980)*100,
                                  ecount(new90), ecount(new90)/ecount(g1985)*100,
                                  ecount(lost90), ecount(lost90)/ecount(g1985)*100,
                                  ecount(new95), ecount(new95)/ecount(g1990)*100,
                                  ecount(lost95), ecount(lost95)/ecount(g1990)*100,
                                  ecount(new00), ecount(new00)/ecount(g1995)*100,
                                  ecount(lost00), ecount(lost00)/ecount(g1995)*100,
                                  ecount(new05), ecount(new05)/ecount(g2000)*100,
                                  ecount(lost05), ecount(lost05)/ecount(g2000)*100,
                                  ecount(new10), ecount(new10)/ecount(g2005)*100,
                                  ecount(lost10), ecount(lost10)/ecount(g2005)*100),
                                4, 8))
colnames(newlosttab) <- c("1975", "1980", "1985", "1990", "1995", "2000", "2005", "2010")

xtable(newlosttab, digits=0)

mean(c(ecount(new75)/ecount(g1970)*100, 
     ecount(new80)/ecount(g1975)*100,
     ecount(new85)/ecount(g1980)*100,
     ecount(new90)/ecount(g1985)*100,
     ecount(new95)/ecount(g1990)*100,
     ecount(new00)/ecount(g1995)*100,
     ecount(new05)/ecount(g2000)*100,
     ecount(new10)/ecount(g2005)*100))

mean(c(ecount(lost75)/ecount(g1970)*100,
     ecount(lost80)/ecount(g1975)*100,
     ecount(lost85)/ecount(g1980)*100,
     ecount(lost90)/ecount(g1985)*100,
     ecount(lost95)/ecount(g1990)*100,
     ecount(lost00)/ecount(g1995)*100,
     ecount(lost05)/ecount(g2000)*100,
     ecount(lost10)/ecount(g2005)*100))


### Table A3. New Embassies Received by Country (1995-2005)

idc00 <- data.frame(sort(degree(new00, mode = "in")), stringsAsFactors = F)
df00 <- data.frame(cbind(row.names(idc00), as.numeric(as.character(idc00[[1]]))), stringsAsFactors = F)
colnames(df00) <- c("stateabb", "embs")
df00$embs <-as.numeric(as.character(df00$embs))
sapply(df00, mode)
df00 <- df00[df00$embs>0,]
df00 <- df00[order(df00$embs, decreasing = T),]

idc05 <- data.frame(sort(degree(new05, mode = "in")), stringsAsFactors = F)
df05 <- data.frame(cbind(row.names(idc05), as.numeric(as.character(idc05[[1]]))), stringsAsFactors = F)
colnames(df05) <- c("stateabb", "embs")
df05$embs <-as.numeric(as.character(df05$embs))
sapply(df05, mode)
df05 <- df05[df05$embs>0,]
df05 <- df05[order(df05$embs, decreasing = T),]

df00
df05


### Network Centrality

### Table A4. Degree Distribution

idc1970 <- degree(g1970, mode ="in")
odc1970 <- degree (g1970, mode = "out")
summary(idc1970)
summary(odc1970)

idc1975 <- degree(g1975, mode ="in")
odc1975 <- degree (g1975, mode = "out")
summary(idc1975)
summary(odc1975)

idc1980 <- degree(g1980, mode ="in")
odc1980 <- degree (g1980, mode = "out")
summary(idc1980)
summary(odc1980)

idc1985 <- degree(g1985, mode ="in")
odc1985 <- degree (g1985, mode = "out")
summary(idc1985)
summary(odc1985)

idc1990 <- degree(g1990, mode ="in")
odc1990 <- degree (g1990, mode = "out")
summary(idc1990)
summary(odc1990)

idc1995 <- degree(g1995, mode ="in")
odc1995 <- degree (g1995, mode = "out")
summary(idc1995)
summary(odc1995)

idc2000 <- degree(g2000, mode ="in")
odc2000 <- degree (g2000, mode = "out")
summary(idc2000)
summary(odc2000)

idc2005 <- degree(g2005, mode ="in")
odc2005 <- degree (g2005, mode = "out")
summary(idc2005)
summary(odc2005)

idc2010 <- degree(g2010, mode ="in")
odc2010 <- degree (g2010, mode = "out")
summary(idc2010)
summary(odc2010)

sort(idc1970)
sort(idc1975)
sort(idc1980)
sort(idc1985)
sort(idc1990)
sort(idc1995)
sort(idc2000)
sort(idc2005)
sort(idc2010)


### Table A5. Power Law Tests (2005)

power.law.fit(idc2005, implementation = "plfit")
power.law.fit(odc2005, implementation = "plfit")


### Table A6. Distributions of Normalized Centrality Measures (2005)

# in- and out-degree
idc <- degree(g2005, mode ="in")
odc <- degree (g2005, mode = "out")

# normalized degree centrality
ndc <- degree (g2005, mode="in", normalized=T)
summary(ndc)

# closeness centrality
cc <- closeness(g2005, mode="all", normalized=T)
summary(cc)

# betweenness centrality (needs to be logged!)
bc <- betweenness(g2005, directed=T, normalized=T)
summary(bc)

# eigenvector centrality
ec <- evcent(g2005, directed=T)
ec <- ec$vector
summary(ec)


### Table A7. Correlations Between Normalized Centrality Measures (2005)

cor(ndc, cc)
cor(ndc, bc)
cor(ndc, ec)
cor(cc, bc)
cor(cc, ec)
cor(bc, ec)


### Figure A1. Degree and Eigenvector Centrality Rankings (2005)

c_rank <- data.frame(ndc, ec, idc)
c_rank$stateabb <- row.names(c_rank)

cor(c_rank$ndc, c_rank$ec)
c_rank$dc_rank <- rank(-c_rank$ndc, ties.method = "min")
c_rank$ec_rank <- rank(-c_rank$ec, ties.method = "min")
c_rank$idc_rank <- rank(-c_rank$idc, ties.method = "min")
cor.test(c_rank$dc_rank, c_rank$ec_rank, method="kendall")

#pdf(file=file.path(SAVEDIR, "DuqueFigA1.pdf"), width=8, height=8)
ggplot(data=c_rank) + geom_text(aes(x=dc_rank,y=ec_rank,label=stateabb)) +
  xlab("Degree Centrality Ranking") + ylab("Eigenvector Centrality Ranking") +
  theme_bw() + theme(axis.title=element_text(size=16), axis.text=element_text(size=14))
#dev.off()


### Figure A2. Hub and Authority Scores (2005)

# Calculate hub and authority scores
hub <- hub.score(g2005)
authority <- authority.score(g2005)

# Check hub and authority scores
summary(hub$vector)
summary(authority$vector)

# Correlations between these scores and/or other centrality measures

ndc <- degree (g2005, mode="in", normalized=T)
nodc <- degree(g2005, mode="out", normalized=T)

cor(hub$vector, authority$vector)
cor(hub$vector, nodc)
cor(authority$vector, ndc)
cor(authority$vector, ec)

#pdf(file=file.path(SAVEDIR, "DuqueFigA2.pdf"), width=10, height=3.75)
par(mfrow=c(1,3))

plot(hub$vector, authority$vector, 
     xlab='Hub score', 
     ylab='Authority score', 
     cex.lab=1.45,
     main='Hub and Authority Scores \nof Diplomatic Network')

plot(hub$vector, nodc, 
     #col=adjustcolor("blue", alpha=0.2), 
     xlab='Hub score', 
     ylab='Out-degree',
     cex.lab=1.45,
     main='Hub Score and Out-Degree \nof Diplomatic Network')

plot(authority$vector, ndc, 
     #col=adjustcolor("red", alpha=0.2), 
     xlab='Authority score', 
     ylab='In-degree', 
     cex.lab=1.45,
     main='Authority Score and In-Degree \nof Diplomatic Network')

#dev.off()


#### Figure 2. Relationship between Material Capability and Degree Centrality (2005)

c_rank <- data.frame(ndc, ec, idc)
c_rank$stateabb <- row.names(c_rank)

vcov <- read.csv("vcov.csv")
c_rank <- join(c_rank, vcov[vcov$year==2005,], by="stateabb")

#pdf(file=file.path(SAVEDIR, "DuqueFig2.pdf"), width=6.5, height=5.2)
ggplot(data=c_rank) +
  stat_smooth(aes(x=log_cinc, y=log(1+ndc)), method="lm", se=FALSE, size = 0.5, colour = "black") +
  geom_text(aes(x=log_cinc,y=log(1+ndc),label=stateabb), size=2.8) +
  xlab("CINC Score") + ylab("Degree Centrality") +
  theme_bw() + 
  theme(axis.title=element_text(size=8.5), axis.text=element_text(size=7))
#dev.off()
